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ABSTRACT 



GRB 050904 at redshift z = 6.29, discovered and observed by Swift and with 
spectroscopic redshift from the Subaru telescope, is the first gamma-ray burst 
to be identified from beyond the epoch of reionization. Since the progenitors of 
long gamma-ray bursts have been identified as massive stars, this event offers 
a unique opportunity to investigate star formation environments at this epoch. 
Apart from its record redshift, the burst is remarkable in two respects: first, 
it exhibits fast-evolving X-ray and optical flares that peak simultaneously at 
t ~ 470 s in the observer frame, and may thus originate in the same emission 
region; and second, its afterglow exhibits an accelerated decay in the near-infrared 
(NIR) from t ~ 10 4 s to t « 3 x 10 4 s after the burst, coincident with repeated 
and energetic X-ray flaring activity. We make a complete analysis of available 
X-ray, NIR, and radio observations, utilizing afterglow models that incorporate 
a range of physical effects not previously considered for this or any other GRB 
afterglow, and quantifying our model uncertainties in detail via Markov Chain 
Monte Carlo analysis. In the process, we explore the possibility that the early 
optical and X-ray flare is due to synchrotron and inverse Compton emission from 
the reverse shock regions of the outflow. We suggest that the period of accelerated 
decay in the NIR may be due to suppression of synchrotron radiation by inverse 
Compton interaction of X-ray flare photons with electrons in the forward shock; 
a subsequent interval of slow decay would then be due to a progressive decline in 
this suppression. The range of acceptable models demonstrates that the kinetic 
energy and circumburst density of GRB 050904 are well above the typical values 
found for low-redshift GRBs. 

Subject headings: gamma rays: bursts: individual: GRB050904 - cosmology: 
miscellaneous 
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INTRODUCTION 



One of the most exciting resu lts from the first year of operations of NASA's Swift 



( ICusumano et al 



satellite mission (Gehrels et al.ll2004T) has been the discovery and observation of GRB 050904 



20061 ) at redshift z = 6.29 (IKawai et al.l 120061 ) . This burst was initially 



detected by the Swift BAT at 01:51:44 UT on September 4, 2005, and was quickly followed up 
by pointed observations with the X-ray telescope (XRT) and UV/optical telescope on Swift, 
and by numerous ground-based facilities. Early afterglow photometry provid ed the first 



indications for a very high redshift for this event, z > 5.3 (IHaisfip et al. 



20051 ) . prompting 



a global observing campaig n that culminated in the spectroscopic observations by Subaru 
that provided the redshift ( Kawai et al. 2006h and enabled t he first investigation of the 
reionization epoch via GRB afterglow light (ITotani et al.ll2006l ). 



Within a day of the burst, the dis covery of an associated prompt I ~ 14 mag optical/NIR 
flash was reported (jBoer et al.l 120061 ). The timing of this flare, which peaks at t ~ 470 s 
after the burst, is coin cident with X-ray (XRT) and 7-ray (BAT) flares observed by Swift 
( ICusumano et al.l 120071) . This is the second such bright optical flash observed so far, after 
GRB 990123 (lAkerlof et al.l Il999l ) , and may even exceed that event in optical luminosity 
(kann et al.ll2007h . 



Subsequent op tical and NIR observa t ions o f the fading afterglow were reported by 
Haislip et al.l (120061 ) and iTagliaferri et al.l (120051 ). and provide evidence for a slow decay 
phase during the first day after the burst, and a jet break at t m 3 days. A campaign of radio 
observations at the VLA yielded multiple detections of the afterglow at 8 GHz, which exhib- 
ited a s low evolution con sistent with a high circumburst density and extremely high kinetic 
energy (jFrail et al.ll2006l ). The X-ray lightcurve observed by Swift exhibits numerous inter- 
esting features, inclu ding the early flare at t ~ 470 s, vigorous X-ray flaring activity, and 
a possible jet break (ICusumano et al.l 120071 ). Finally, observation s with the Hubble Space 
Telescope and the Spitzer Space Telescope by Berger et al. ( 2006 )^1 have yielded late-time 
detections of the afterglow and host galaxy. 



Studies of GRBs at low-redshift, l<,z<,3 (IPanaitescu &: Kumarll2001al : lYost et al.ll2003l ). 
have established that the circumburst densities for these events range between 0.06 and 
30 cm -3 , with their kinetic energies having a relatively narrow dis tribution with a peak 
around 5 x 10 50 ergs (IPanaitescu fc Kumarl 120021 : iBerger et al.ll2003l ). It is therefore inter- 
esting to investigate whether the quantities for the high-redshift GRBs follow these results 



1 Bergen (|2007f ) have re-observed the position with Hubble Space Telescope on July 22 UT, 2006 and then 
provide the upper limit on the host galaxy. In the modeling, we have used the updated data point. 
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or not. 

Several questions seem particularly pertinent. What density range is found in the en- 
vironments of GRBs at high redshift? How do their kinetic energies compare to those of 
low-redshift events? Are other properties of high-redshift GRBs - their beaming angles 
and shock microphysical parameters - the same as for low-redshift GRBs? The answers to 
these questions can potentially cast light on outstanding mysteries of the GRBs themselves, 
and reveal important aspects of the early Universe. These interesting questions provide the 
motivation for our detailed investigation of the properties of GRB 050904. 

In this paper, we attempt a complete model of the full set of X-ray, optical/NIR, and 
radio afterglow observations of GRB 050904. Our derivation of physic al parameters from 



afterglow observations is carried out in the context of the fireball model (|Meszarosll2006l . and 
references therein). 

We make explicit consideration of two scenarios for the origin of the X-ray and optical 
flares at t fa 470 s. Our first scenario (A) attributes the flares to internal shocks or engine 
activity, and excludes the flares from a fterglow fits, wit h only later observations considered. 



This is consistent with the approach of I Wei et al.l (120061 ). who argued on the basis of the fast 



decay of the flares that they could not be due to reverse shock emission. In our alternate 
scenario (B), however, we attribute the flares to emission from the reverse shock. In order to 
accommodate their fast decay, we use a starting time for the asymptotic Blandford-McKee 
solution which is near the start of the flare, later than the burst trigger time, which serves 
to flatten the post-flare decay. In this scenario, the optical flare comes from synchrotron 
radiation in the reverse shock, and the X-ray and 7-ray flares are produced by synchrotron 
self-Compton scattering (SSC) of photons in the reverse shock region. A cosmology where 
H = 71 km s _1 Mpc -1 , Q\ = 0.73, and fl m = 0.27 is assumed in calculating the luminosity 
distance D L . For GRB 050904 at z = 6.29, the luminosity distance is D L = 1.93 x 10 29 cm. 

Our model and its supporting analytical formulae are presented in §2j with several 
derivations reserved for the appendices. Our numerical simulation procedure and the data 
set used in our fits are described in §31 while in §H we analyze the results, including the 
J-band light curve ( §4.ip . the radio light curve ( §4.2p . and our ultimate energy and density 
constraints ( §4.31) . A discussion of the results and our conclusions are presented in £j5] and 
§61 respectively. 



-4- 



OBSERVATIONS AND THEORETICAL FRAMEWORK 



2.1. Burst and Afterglow Observations 

The Swift BAT observations of the prompt emission give a 7-ray duration of Tg = 
225 ± 10 s, a spectru m with power-law phot on index T = 1.34 ± 0.06, and a fluence of 



5.4 x 10 6 erg cm 2 (jCusumano et al.ll2006l ). Given the burst redshift of z = 6.29 ± 0.01, 

10 54 ergs. XRT observations began 



the isotropic-equivalent gamma-ray energy is E, 



7, ISO 



161s after the burst and continued for 10 days after the burst trigger, overlapping with BAT 
observations for about 300 s before the high energy emission faded below the BAT threshold 
(jCusumano et al.ll2006l ). 



Thanks to the Swift prompt alert, the early afterglow was also observed promptly by 
the TAROT robotic telescope. The TAROT observations started at 86 s after trigger and 
lasted for more than 1500 s; by making a spectrophotometric calibration of the field, Boer et 
al. (2006) are able to present their data as flux densities at 9500 A, which we use and shall 
refer to as the TAROT /-band o bservations. Other l arger ground-based tele scopes started 



imaging the field 3 hours later (IHaislip et al.l 120061 ; iTagliaferri et al.l 120051 ) . We present 



our compilation of the observational data in the X-ray (XRT), hard X-ray (BAT), and 
optical/NIR in Fig. [IJ Here and throughout the paper we convert X-ray measurements 
to flux density measurements at the frequencies of 5 keV (XRT) or 50 keV (BAT); these 
energies correspond to observing frequencies of 1.2 x 10 18 Hz (XRT) and 1.2 x 10 19 Hz (BAT), 
respectively. The conversion factors from photon counts per second (cps) to flux are 3.31 
/iJy cps" 1 for XRT PC data, 1.82 ^Jy cps" 1 for XRT WT data, 154.6 /iJy cps" 1 and 86.2 



/iJy cps 1 for BAT masktag-lc data, respectively (see ICusumano et al.l 120071 for the details 
of XRT PC, WT and BAT data). 

As seen in Fig. [TJ the flux and spectral evolution of the burst and afterglow divide 
the lightcurve into six distinct segments, determined by inspection and motivated by the 
physical model put forth herein. These segments are: (A) t < 350 secon ds: In the X-ray 
band , the flux decays as F v oc t~ a with an index of a = 2.07 ± 0.03 (jCusumano et al. 



20061 ). We follow the conventional definition for the flux F v oc t a v 13 . In the J-band 



observation, there are two observational data points, the earlier of which is only an upper 
limit, apparently indicating an increasing tendency of the flux with time. (B) 320 < t < 600 
seconds: Flares are observed in both the / and X-ray bands, and also in the BAT energy 
range. They peak around 470 seconds after the burst trigger time. The spectral index 
evolves from 0.50 ± 0.07 to 0.88 ± 0.12 over this time interval. (C) 600 < t < 6.3 x 10 3 
seconds: A power-law decay is shown in the X-ray band. During the same time interval, 
there is no optical/NIR detection, except for two upper limit flux values at 9500 A. (D) 
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6.3 x 10 3 < t < 4.3 x 10 4 seconds: Many irregular fluctuations are observed in the X-ray 
band . The J band shows a decay which can be described with a = 1.36 +0 "' 

4.3 x 10 4 seconds) < t < 2.6 days 



q Qg (lHaislip et al. 



2006). (E) 0.5 days 



There is no effective XRT 
observation within this period, and no further fluctuations are detec ted. The flux in th e J 
band is flatten ed a little bit with a t emporal index of a = 0.82to!os (lHaislip et al.ll2006l ) or 
a = 0. 72^20 (ITagliaferri et al.ll2005l ). (F) t > 2.6 days: the flux decay index in J band is 
around a = 2.4 ± 0.4. Only one data point in the X-ray band is available. The J-band data 
shows a sharp break, which is thought to be the jet break. 



2.2. Afterglow Modeling in the Swift Era 

Thanks to the rapid and precise alerts generated by Swift, and its extensive multiwave- 
length follow-up campaigns, the afterglow data collected during the Swift era has resulted in 
some necessary modifications to the standard afterglow model. 

Data from Swift have provided the greatest advance over earlier datasets in the X-ray 
band, as Swift responds to the burst orders of magnitudes faster than previous generations 
of satellites and can often track the X-ray afterglow for up to 10 days. Many new features 
of the X-ray afterglow have thus emerged, leading to the identification of a canonical X-ray 
afterglow behavior. In addition to the prompt emission phase, this involves at least five 
components of the X-ray afterglow (Nousek et al. 2006; Zhang et al. 2006a), which are: 
(1) A steep decay phase, often interpreted as the tail of the prompt emission, and thus, 
part of the GRB internal shock; (2) A shallow decay phase of uncertain origin, with several 
theoretical models proposed, including energy injection, jet inhomogeneities, or varying shock 
microphysical parameters; (3) A normal decay phase, familiar from pre-Swift observations; 
(4) A post jet break phase; and (5) The X-ray flares, which are superposed on the various 
power-law segments of the afterglow's decay, and are fast-evolving in the sense that the 
flare's rise and decay timescales 5t are much smaller than the time since the burst t, that 
is, dt/t <C 1. The current interpretation ascribes the X-ray flares to the same cause a s 



the prompt gamma-ray emission - energy dissipation in internal shocks (jZhang et al.ll2006l ). 



Since the X-ray flares are thus a manifestation of central engine activity, it is necessary 



to exclude them from afterglow model fits (e.g., iFalcone et al.ll2006l ; IChincarini et al.l 12007 



Nava et all 120071 ). 



The accumulation of Swift afterglow observations has also raised questions about the 
collimated or "jet" interpretation of afterglow lightcurves. Traditionally, jets have been 
invoked to explain a late time (£>1.0day), broadband and achromatic steepening in the 
afterglow decay by an increment of +1 in the power-law index, from a ~ 1.2 to a ~ 2.2. 
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Specifically, the late-time decay index is predicted on robust grounds to be a = p, where p 
is the power-law index of the energy distribution of the synchrotron-emitting electrons that 
generate the afterglow light. 

Current challenges to the jet break picture are two-fold. Fir st, there are few jet break s 
seen in the X-ray light curves of afterglows observed by Swift (IBurrows fc Racusinl 120071 ). 
As of 7 October 2006, the XRT had observed 145 long GRBs and 16 short GRBs, with 
the X-ray afterglows of almost all long bursts followed for days to weeks with Swift, in 
campaigns that typically last until 10 days after the burst. Unexpectedly, among these 
bursts, only 4 long burst s and 2 short bursts show the jet break signature (See Table 1, 
Burrows fc Racusinl 120071 ) . Second, where X-ray breaks of the appropriate steepness are 
seen, they are often observed to be chromatic, exhibiting a different evolution in UVOT (or 
ground-based optical) observations. Thus the jet break picture is being reevaluated, with 
new possibilities, such as en ergy injection or t he evolution of microphysical shock parameters, 
under active consideration (jPanaitescul 120051 ) . 



Despite these results, we believe that the jet break picture remains the better inter- 
pretation for most bursts, even in the Swift era. First, the energetics of GRBs are ex- 
tremely difficult to reproduce from stellar-mass progenitors without beaming corrections of 
some sort. Second, the successes of the standard jet break picture in modeling bursts from 



before Swift are too numerous and significan t to be discounted (e.g.. lHarrison et al.l 11999 
Panaitescu &: Kumarll2001bl ; lYost et al.ll2003l ). Third, several candidate jet breaks have been 
ident ified in the Swift era, and when these are seen, they exhibi t the expected properties 



(.e.g. 



Cenko et all 120071 : iDai et al.ll2007l : IBurrows fc Racusinl 120071 ) . 



Finally, several arguments support the presence of a jet break in the evolution of the af- 
terglow of interest to us here, GRB 050 904. The break that is o bserved in the NIR lightcurve 
of the afterglow at t — 2.6 ± 1.0 day (ITagliaferri et al.l 120051 ) is achromatic, being seen in 
multiple NIR bands, and is followed by a steep post-break decay with power-law index 
a = 2.4 ± 0.4, consiste nt with the closure relation expected from the sta ndard jet model 
(Ta gliaferri et al.l 120051 ). Separately, combi ning the Ghirla nda relation (IGhirlanda et al. 



20041 ). and expression for the jet break time (jSari et al.lll999l ). a jet break is expected in the 
range between 0.64 and 3.0 days. The peak energy of GRB 050904 is E s ™ ak > 150(1 + z) — 
1095 keV and its isotropic-e quivalent gamma-ray energy lies between 6.6 x 10 53 and 3.2 x 10 54 
erg (jCusumano et al.l 120061 ) . so we take E^ k ~ 1100 keV and -E^,,;^,,^ ~ 100- lu addition, 
the circumburst density lies between 40 < n < 1000 cm -3 (IKawai et al.l 120061 ) . Setting 
n rsj 100 cm -3 and the radiation efficiency in the prompt emissi on phase V y = 0. 2, one can 



obtain the jet break ranges of 0.64 and 3.0 days from Eqn. 3 of ISato et al.l (120071 ). which is 
consistent with the observed break time. 
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2.3. Two Different Scenarios 



2.3.1. (A) Forward Shock Fit Only 



It has been argued by IWei et al.l (120061 ) that the flares at t ~ 470 s are due to internal 
shocks, rather than a reverse shock, based on an apparent steep decay of the X-ray light curve 
right afte r the peak, with a tempora l index of a± ~ 8.8 when referenced to the burst trigger 



time (see IKobayashi fc Zhang) 120071 for a discussion on the onset of GRB afterglow). This 
is because the fastest decay must not be steeper than that from the high latitude emission, 
whose temporal index is a<i = (3 + 2, in this case, a ~ 3. This leads them to favor internal 
shocks for both the optical and X-ray flares, or a combination of reverse shock for the optical 
and internal shocks for the X-ray. Consistent with this argument, we investigate a model of 
the forward-shock emission only, that is, we consider only data points in regions (C)-(F) of 
Fig. 1. 



2.3.2. (B) Forward Shock and Flares Fit 

Separately, we consider a scenario in which the flares peaking at t ~ 470 seconds in 
the optical and X-ray are due to the reverse shock. This is motivated by the following 
argument. First, from the observational point of view, the flares of GRB 050904 have a great 
simil arity with the beha vior of GRB 990123, which is thought to be due to the reverse shock 



(e^. lSari fc Piranlll999h . Second, the Blandford-McKee (BM) solution flBlandford k McKee 
19761 ) describing the afterglow is parameterized with a time origin to given by the trigger 
time. However, the BM solution is an asymptotic description of the afterglow evolution, 
which is valid for times substantially longer than any of the timescales associated with the 
prompt emission. Clearly, however, there is a transition from the initial prompt phase, when 
the bulk Lorentz factor is relatively constant, to the steady deceleration phase when the 
asymptotic BM solution applies. Numerical simulations have so far not provided specific 
guidance on the most appropriate value of the reference time for BM (afterglow) evolution. 

Naturally, however, the steepness of the light curve decay, as parameterized by F v ~ 
{t — t )~ a , depends on the reference time t that is adopted. This has been discussed most 
recently in the context of the X-ray flares s een by Swift, whi ch exhibit steep decays and 



so are widely attributed to internal shocks (ILiang et al.ll2006l ). Here we are dealing with 



a somewhat different situation, in that the X-ray light curve exhibits an initial flare and 
a subsequent decay, which we propose to attribute to a reverse shock. In the absence of 
guidance from numerical simulations, we adopt a purely phenomenological approach, based 
on the constraint that, whether for internal shocks or for a reverse external shock, the 
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temporal decay cannot be steeper than (t — t' ) a , where t' Q is the reference time which 
fits the high latitude decay a = 2 + /3, where j3 is the spectral inde x. For late centra l 



engine internal shocks, t is found to be near the onset of the last spike ( Liang et al.ll2006l ). 
Here, in the same spirit, since the initial X-ray flare is assumed to be due to the reverse 
shock, after which the BM solution is asymptotically reached, the new reference time t is 
determined from the constraint that the decay slope be equal to that expected from the high 
latitude mechanism. Having set a new reference time t , all the time- related quantities for 
the afterglow evolution, such as u m , u c and so on, are now referred to (t — 1' ). We emphasize, 
however, that although the reference time is shifted from the burst trigger time to a new 
time point, the physical deceleration timescale remains the usual one. For GRB 050904, the 
best-fitting reference time satisfying the high latitude condition is t = (0.86 ± 0.01)^0 
where the (usual) deceleration timescale is 468 s (and the latter is measured relative to the 
burst trigger time). 

The deceleration time is a critical point in the afterglow evolution. After that, both the 
reverse shock and the forward shock evolve into the asymptotic regime of the BM solution. 
In model (B), the peak time should represent the deceleration time of the shock. In the 
discussion below all the formulae after the deceleration time are taken relative to new refer- 
ence time t' (we have inserted the exact value of t into the equations below). A difference 
between models (A) and (B) lies in that t' = for model (A) and t > for model (B). The 
determination of t' for model (B) will be described in §3.11 



2.4. Synchrotron and Inverse Compton Afterglow Model 



The flares in the optical (jBoer et al.ll2006l ) and X-ray (jCusumano et al.ll2006l ) are ob- 
served to peak simultaneously at t ~ 470 seconds after the burst. The X-ray data has better 
time resolution, with an estimated peak time of t pea k — 468.0 ± 2.0 s. 

We consider fits for two generic models, as described above. In model (A) the flares 
are ignored, as they may come from other regions (possibly internal shocks) rather than the 
reverse shock region. In model (B), however, we consider the flares to be due to the reverse 
shock. In this case, since the peak of the flare is c learly separated fr om the GRB itself, 
the reverse shock must be in the thin-shell regime (jZhang et al.l 120031 ) . Spectral analysis 
(jBoer et al.l 120061 ) shows that the flares in the X-ray and optical bands must be due to 
different mechanisms. Here, we assume that the optical flare is due to synchrotron radiation 
in reverse shock region, and that the X- ray flare is produced by synchrotron-self- Compton 



(SSC)in the reverse shock region (e.g., IWang et al.l l200ll ; iKobayashi et al.l 120071 ). Other 



observed X-ray components at late times, except for the flares, are assumed to be produced 
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by synchrotron radiation from the forward shock. We use the following reference values for 
the main parameters involved: the isotropic-equivalent kinetic energy E 52 = 1, the external 
circumstellar density n = 1 cm -3 , and the magnetic field ratio of the reverse and forward 
shocks R B = eB,r/^B,f- In addition, we have assumed equality of the energy equipartition 
parameter for electrons in the forward and reverse shock regions in model (B). 



A transition from the temporal index a = 0. 82 ± 0.08 to a = 2.4 ± 0.4, interpreted as 
the jet break, is observed at £j et = 2.6 ± 1.0 days (ITagliaferri et al.l 120051 ) . and the sharpness 



of this break makes the burst not likely to have occurred in a stellar wi nd-type density 
environment (IKumar &: Panaitescul l2000l ; iGou et al.ll200ll ; iFrail et al.ll2006l ). Therefore, we 
will focus on afterglow evolution in a uniform density environment throughout the paper. 



2.4-1- Forward Shock Synchrotron Formulae 



Usually the synchrotron emission spectrum of the forward shock is described by a bro- 



ken power-law with three critical quantities z/„ 



and -F^ ax which are the synchrotron 



frequencies for the electro n energies at the injection, cooling and peak flux emission, re- 
spectively (jSari et al.lll998l ). In modeling the radio emission, it is necessary to also take into 
account the self-absorption frequency, v sa . We present the formulae including self-absorption 
in Appendix |Al because there ar e different forms dep ending on various afterglow regimes, 
six in all, for the spectra given by lGranot fc Saril (120021 ) . The main quantities of interest are: 



2.65 x 10 23 Hz 



P 



p — 1 



1/2 2 < . _ . 
€ Bj e e \ l l 



rpsyn 
u,max,f 



1.58 x 10 15 Hz E^ezfnQ 1 [t-t 



,x -3/2 fl + Z 

V 

,\-l/2 



7.29 



1/2 



X 



[l + yy 



l + z 
7.29 



-1/2 



7.63 x 10" 5 Jy e^ 2 f E 52 n /2 



Dr(z) 



1.9 x 10 29 



-2 



where the convention Q = 10 X Q X is used, kinetic energy E = lO 52 -^^ erg, and density 
n = 1 n cm -3 . e B is the magnetic field equipartition parameter, and e e is the electron 
equipartition parameter. The subscript "/" denotes the forward shock, and "r" denotes the 
reverse shock. The parameter Y re fers to the first-ord er inverse Compton effect and is defined 
as Y = (—1 + a/1 + r]e e /eBj)/2 (jSari fc Esinl 120011 ). Here r/, the fraction of the electron 
energy that is radiated away, expresses the magnitude of the radiative correction: i] = 1 for 
fast cooling and rj = ('Jc/lm) 2 ~ p for slow cooling, where 7 C is the cooling Lorentz factor and 
7 m is the typical Lorentz factor for the electrons. For z = 6.3, the luminosity distance in a 
concordance cosmology is ~ 1.9 x 10 29 cm. After the deceleration time, relative to the 
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new reference time, the afterglow will evolve asymptotically as the BM solution. Thus, the 



evolution relation for each quantity is u c oc (t—t 



-1/2 



The quoted typ ical values for u„ 



and 



Ff^ a f are consistent with the ones in 



Granot fc Saril (120021 ). but the values for v m an d v r are 2 times smaller, and the typical 



value for F^ ax j is 4 times smaller than that in ISari et al.1 (119981 ). which are fitted for the 
numerical simulation results. 



2.4-2. Reverse Shock Synchrotron Formulae 



In the reverse shock at the deceleration time, the typical quantities are (IKobayashi et al. 



20071): 



— RBRMF v , max j (2) 

where R B = (e B)T / e BJ ) l l\ R x = (1 + Y)/(l + X), X = (-1 + ^1 + 7^/^0/2 is the 
Compton parameter in the reverse shock region, and the factor Rm = Trf/To, where Td is the 
Lorentz factor at d eceleration time. For the thin-shell case, « r , so we have Ym ~ T 



(IZhang et all 120031 ). 



After the deceleration time, since the shock is in the thin-shell case, each quantity 
evolves as v c oc (t — t' ) ac where a c = (15 g + 24)/(14g + 7), v m oc (t — t' ) am w here a m - 



-(15 g + 24)/(Ug + 7), and F vm oc (t-t' ) a f where a f = -(llg + 12)/[7(1 + 2g)} (IZou et al. 



20051 ) where (7 is the evolution index for r oc R 9 where V is the Lorentz factor of the 



afterglow, and R is the radius of the afterglow. 



2-4-3. Inverse Compton Effects, Jet Break and Non-relativistic Case 



Both models include inverse Compton effects. The formulae f or the norm a l case v a < 
v m < v c and v a < v m < v c in the forward shock case are listed in ISari fc Esinl (120011 ). We 
have derived the formulae for additional cases where the self-absorption frequency is above 
the typical frequency v m or cooling frequency v c in Appendix [B] For inverse Compton effects 
in the reverse shock, the self-absorption frequency has a similar form as that of the forward 
shock, the only difference being that the forward shock-related quantities are replaced by 
the corresponding reverse-shock quantities. 

Before the bulk Lorentz factor drops below the inverse of the opening half-angle, 
each of the typical quantities follows the evolution given above (see sectio ns 12.4. II and |2A2p . 
After T < 6~ l , those typical quantities follow a different evolution given by lSari et al.1 (119991 ). 
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We also consider the non-relativistic (NR) evolution of the afterglow in its end-state. 
The time for the afterglow to enter into the NR stage is calculated by the condition that the 
Lorentz factor of the shell T = 2. After the afterglow evolves into the NR stage, it s dynamics 



are de scribed by the self-similar Sedov-von Neumann- Taylor solution, for which iFrail et al. 
( 120001 ) have given a detailed description. 



2.4-4- Host Galaxy Extinction and Lyman-a Damping Absorption 

We have considered the host galaxy extinction in the optical band in the rest frame of 
the host galaxy, using a Milky Way extinction curve. We note tha t the particular ty pe of 



extinction curve used is not expected to make a difference in this case; iKann et al.l (120071 ) have 
tested the application of SMC, LMC, and Milky Way extinction curves to the composite J- 
band data, and all three models suggest minimal extinction. In our fitting, we treat the visual 
extinction Ay in the host galaxy as a free parameter. In addition, because the extinction 
affects the spectral shape, we have considered the spectral index correction due to the host 
galaxy extinction in our fitting. 

Besides the normal extinction by the host galaxy, the emission close to the wave- 
length of Lyman-a at z = 6.29, including z band emission (effective wavelength 0.91 fim 
and band w i dth is 0.13 fim), has undergone neutral hydrogen absorption in the IGM. 



Totani et al.l (120061 ) fit the spectrum of GRB 050904 and find its best-fit column density 
to be log JVjj/(cm' 2 ) = 21.62. To find the expected Lyman-a absorption in the z band, 
we have convolved the Lyman-a absorption profile with the filter transmission curve of the 
z-band, and obtain the absorption coefficient A = 0.77 meaning 33% loss of z-band flux. 



Also we notice that from Figure 6 in lTotani et al.l (120061 ). the absorption around wave- 



length 9500A is negligible, so we take the absorption efficient A = 1 for the 9500 A data. 

3. Fitting Data and Procedure 

3.1. Observational Data and Constraints 

Our model fits for the two different scenarios cover a range of bands from the radio, 
through the IR/optical, to X-ray and BAT energies. The list of the observational data used 
in our global fitting is in Table HJ 

In model (A), the new reference time is the trigger time, so to = 0. In model (B) with the 
flares included, we have to determine the new reference time first. We will now describe our 
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procedure for the early time data used as input to the fits of model (B). Relative to the trigger 
time, the temporal index of the fast decay is a ~ 8.8. The fast decay is considered to be 
due to the high latitude emission of the reverse shock. In order to obtain the new reference 
time, we assume a simple power-law model for the high latitude flux f v = F u ( tobs ~ t{, )~ a , 
where F v denotes the peak flux at the peak time t p , t is the starting time at which the 
asymptotic Blandford & McKee solution applies, and a is the temporal index. The high 
latitude emission is assumed to decay with an index ah = 2 + (3. For the X-ray observations 
right after the peak, t > 470 seconds after the burst, the observational spectral index is 
around (3 = 0.88 ± 0.12, so we take the temporal index in the fitting as 2.88 ± 0.12, and 
selecting two observed data points, we can find the new reference time to be (0.86±0.01)t,i ec . 
In our fitting, we have set the new reference time at 0.86td ec . 

We notice that right after the fast decay starting 470 seconds after the burst there is a 
plateau, extending from t ~ 600 seconds to t w 2000 seconds (see Fig. [2]). A close look at 
the data around 1000 seconds shows an obvious flare, so we ignore the data after that. Thus, 
we only keep the X-ray data between 470 seconds and 1000 seconds. For this part of the 
data, there are three different contributions: (1) Forward shock; (2) High-latitude emission 
of the flare; (3) High-latitude emission of the prompt emission. It should be pointed that 
the high latitude emission of the prompt emissio n must be described by a broken power-law 
due to the spectral evolution of the burst itself (jCusumano et al.ll2006l ); the break time in 
the rest frame is at t m 350 s in the observer frame. Before the break time, the temporal 
index is a ~ 2, and afterwards, it is a ~ 3. We subtract the high latitude emission from the 
observed flux to obtain the "pure" forward shock emission. To reduce the uncertainty of the 
data points, we then find the mean flux by averaging adjacent data points in groups of five, 
and use the averaged value for our fitting, finally providing the fourth data point in Table [TJ 

We also apply a similar subtraction method to the data around t ~ 2000 seconds 
indicated in Fig. [2J The averaged data point is the 5th row in Table [TJ One difference 
between the subtraction for the fourth and fifth data points is that we only consider the flare 
contribution for the fifth data point (the prompt emission is considered negligible here), 
while both the flare and the high latitude emission of the prompt emission are considered 
for the fourth data point. 

The common constraints considered in both model (A) and model (B) are the following: 
(1) The jet break ti me is fitted to be 3.17 ± 0.22 days (1-er) by combining all available 
NIR and X-ray data. iKann et al.l (120071 ) have extrapolated all the other NIR/optical data 
to the J ba nd and made a com posite light curve in the J band, including the HST data 
observed by iBerger et al.l (120061) ~ 23 days a fter the burst, and obtain a jet break time 



tj e t = 2.63 ± 0.37 days. iTagliaferri et al.l (120051 ) give a value of tj et = 2.6 ± 1.0 days based on 
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multiband fi tting of a smaller data set. (2) The spectral index in the J-band at 1.155 days is 
1.25 ± 0.25 fjTagliaferri et al. 2005 ). (3) The average spec tral index for the early time X-ray 



afterglow from t = 680 s to t = 1600s is = 0.96 ± 0.19 flCusumano et all 12009 ) 



Besides the constraints above, we introduce another constraint for model (B): (4) the 
deceleration time is tdec = 468.0 ± 2.0 s from our estimate of the peak of the X-ray light 
curve. 

It should be mentioned that in model (B) we have summed over both the reverse and 
forward shock flux to fit the observed data (/, X-ray and BAT bands) at the deceleration 
time. 

Thus, considering all the available data listed in Table dj as well as the observed spectral 
index and the jet break time, we have 37 constraints for case (A) and 33 constraints for case 
(B). 

For the fitting and (simultaneous) evaluation of parameter uncertainties, we perform a 
Markov Chain Monte Carlo analysis (MCMC; see Sec. 13.21 for a more detailed description). 
In order for the code to spend most of its time in the regions which have physical meaning, 
we provide penalty conditions during the calculation of the chi-square. The penalty is levied 
by giving an extra boost to the chi-square value when the penalty condition is violated (when 
the condition is satisfied, this extra chi-square is zero). To ensure the smoothness of the fit 
function at the penalty boundary, the penalty function we choose (more or less arbitrarily) 
is A% 2 = [(x — xii m )/(0.01*min(x, xi im ))] 4 , where xn m is the critical value for each parameter. 

We have included some upper limits in our dataset by converting the upper limits into 
synthetic measurements with error bars, as follows. Assuming all the upper limits as 2-cr 
limits (we n otice that the J band data labeled as No. 13 in Table [T] is provided as a 3-cr 
upper l imit, Tagliaferri et al.l2005 ; and the confidence level for the I band data, No. 8, is not 



stated, iBoer et al. hood ), we take the measurement to be half the upper limit, and the error 
bar to be one-quarter the upper limit (half the synthetic measurement). For the radio data, 
which are provided as measurements with error bars even when no detection is realized, we 
adopt these measurements and error bars directly, while plotting 2-cr upper limits in our 
figures. 

For both models (A) and (B), the parameter ranges are restricted as follows: (1) e e < 0.5; 
(2) 6B, r < 0.5; (3) 6bj < 0.5; (4) Electron energy power-law index p > 2.06; (5) 7 m > 2.1 at 
t = 10 7 seconds after the burst. Violations of the parameter ranges incur chi-square penalties 
as discussed above, on a parameter by parameter basis. 

Note that the critical electron energy index value (4) p = 2.06 is found by equating the 
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minimum electron Lorentz facto r for the y > 2 and p = 2 cases. For p > 2, the typical 
Lorentz factor 7 m = j^j^e e r jSari et al.|[l998|). For p = 2, 7 m « ln f4 %j£e e r 0. 



In(4.0xl0 7 ) m e 



Also note that, because we have assumed that the critical radiation mechanism is syn- 
chrotron radiation, the electrons in their shell frame are relativistic, and correspondingly 
7 e > 2 or 7 m > 2 (because at late time the afterglow are in the slow-cooling regime, the 
Lorentz factor of the most electrons are concentrated at 7 m ). We have set the critical value to 
be (5) 7 m = 2.1 rather than 7 m = 2 to eliminate the unrealistic parameter set corresponding 
to 7 m < 2. 

We have not constrained the extinction parameter Ay, since in principle it can be any 
positive number; however, in order to avoid considering negative values of the extinction, 
the code makes use only of the absolute value of this quantity. 



3.2. Parameters and Methodology 



For model (A), we have 7 free parameters, which vary freely, subject only to the penalty 
conditions. The parameters are: the energy index p, the isotropic-equivalent kinetic energy in 
units of 10 52 ergs E 52 , the energy fraction in electrons in the reverse shock and forward shock 
regions e e (note that the electron equipartition parameter in the forward shock is assumed 
to be the same as the one in the reverse shock), the magnetic field equipartition parameter 
in the forward shock €bj, the circumburst density n, the jet opening half-angle 9, and the 
extinction parameter in the host galaxy Ay. In model (B), we introduce two additional 
parameters: the magnetic field equipartition parameter in the reverse shock e# ir , and the 
initial Lorentz factor r . We have assumed that the magnetic equipartition param eter in the 



forward and reverse shock s can be different, whic h is m otivated by the results of iFan et al. 
f l2002h : Izhang et all (jioiii ): iKumar fc Panaitescul (|2003h . 



In order to obtain best-fit parameters and explore the parameter space of the fit func- 
tion, we have tested both a grid search method and a Markov chain Monte Carlo (MCMC) 
method. However, we choose the MCMC after some tests. The grid-based likelihood analysis 
calculates chi-squared values at each grid point of the parameter space, and determines the 
best fit parameters and confidence levels by finding the minimum chi-square value point and 
range of values within a certain "height" above that minimum. The benefit of this method 
is primarily that it is straightforward. Once the parameter ranges and the number of grid 
points are defined, the code is easily implemented. However, the drawback is that it requires 



The typical Lorentz factor for p — 2 case can be obtained in the similar way as the one for p > 2 
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prohibitive amounts of time, especially if there are many free parameters. For example, a 
coarse grid with 7 points per dimension and with 8 parameters requires 2.1 x 10 7 evalua- 
tions, and at 0.2 s per evaluation, the calculation takes r ~ 5 days on a single processor 
machine. Increasing the number of parameters, much less increasing the number of grid 
points, quickly becomes infeasible. By contrast, the MCMC method is very efficient, with 
execution time scaling linearly with the number of parameters, which allows us to perform 
likelihood analyses in a reasonable amount of time. 

Briefly, the MCMC is a method to reproduce, directly, the posterior dis tribution of 



the m odel parameters (for a detailed treatment in an astronomical context, see lVerde et al. 



20031 ). After a limited "burn-in" phase, it should generate a random draw from the posterior 
distribution for most new function evaluations. From this sample, we can then estimate all of 
the quantities of interest for the posterior distributions (the mean, variance, and confidence 
levels). As mentioned above, the MCMC method scales approximately linearly with the 
number of parameters, allowing one to perform a likelihood analysis in a reasonable amount 
of time for a large number of parameters. After an initial burn-in period, and assuming that 
convergence of several chains can be established, all samples can be thought of as coming 
from a stationary distribution. In other words, the chain has no dependence on the starting 
location (although a good choice of starting points and step size can accelerate the chain 
convergence). 

In implementing a MCMC approach, two key and interrelated questions are: (1) At 
what point does the chain converge, that is, how fast does the chain realize the target 
distribution?; and (2) Does the chain provide g ood mixing, that i s, has the chain covered 



all interesting portions of the parameter space? iGelman fc Rubinl (jl992j) suggest a method 



to test t he convergence an d mixing and introduce for this purpose a parameter labeled R 



(see also I Verde et al.ll2003l ). The convergence can be monitored by calculating R for all the 
parameters in two or more chains, and running the simulations until all R values are less 
than 1.2. More conservatively, we may choose to run until all R < 1.1; this is the criterion 
we adopt as our test of convergence. 

As we mentioned above, the MCMC model efficiently explores the parameter space, 
which guarantees that the global minimum will be approached in the long run. By contrast, 
in a grid search method, one has to provide the parameter range beforehand, and one is never 
sure that the minimum chi-square found is the global minimum instead of a local minimum. 
This can lead to very different results. 

Before running the MCMC, it is necessary to initialize the starting point and assign step 
sizes for each parameter. For the starting point, we run a test chain first, then choose the best 
parameter set (evaluated by a minimum chi-square) as the starting point for a formal run. 
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Initially, we set the step size for each parameter to be the l-a range from this same initial 
run; however, several experiments convinced us that a half-sigma step size provided better 
convergence speed. Our MCMC code was implemented in a Matlab environment on a single 
processor machine, and then transplanted to the High Performance Computing (HPC) Linux 
cluster at Pennsylvania State University. We made use of 4 processors, each running one 
chain, and with each chain set to run for 2 million steps at a time. After the chain calculations 
are completed, we merge them and test for convergence. If the chains have not converged 
over the final 1 million steps of the 2 million-step chains, then we take the final parameter 
set as the starting point for another run, execute 2 million additional steps, and test for 
convergence again. The result for our final model, presented here, provided convergence 
after the first run in both cases (Models A and B), with R < 1.1 for all parameters after 
2 million steps. 



3.3. Numerical Results 

For the results presented here, we ran four chains for 2 x 10 6 steps for each of the models 
(A) and (B). Convergence was tested and parameters quantified using the final one million 
steps only. We found that both chains had already converged after the first run. For model 
(A) the R values are 1.008, 1.068, 1.045, 1.022, 1.076, 1.070, 1.004 for the parameters p, 
e e e Bj n, 9, E 52 , A v , respectively. For model (B), the R values are 1.002, 1.004, 1.002, 
1.004, 1.005, 1.005, 1.000, 1.004, 1.000 for the parameters p, e e e B j n, T , 6, e B , r , E 52 , A v , 
respectively. 

The posterior distributions for the parameters for model (A) are displayed in Fig. [3J 
The shaded blue regions delimit the l-a range (68.2%-confidence), and the region included 
within vertical blue lines corresponds to the 90%-confidence interval. We choose the ranges 
of minimum width for these confidence intervals. Separately, we indicate the posterior distri- 
bution of model parameters for models with 7 m < 2.1 (at t = 10 7 s) in green, with a red color 
for its la region. The number of model realizations that were constrained in this sense is 
roughly 1.5% of the total trials. The reduced chi-square for model (A) reaches its minimum 
value 36.2/26 ~ 1.39, and the best fit parameters are p = 2.152, e e = 0.031, e B j = 0.198, 
n = 84.4 cm" 3 , 9 = 0.128, E 52 = 22 A, and A v = 0.034 mag. 

The light curve for the best fitting parameters is shown in Fig. H] (indicated with solid 
line). We have shown the observational data as the background in the grey color, and 
indicated the data actually used for fitting in other bright colors (for detailed description on 
which color stands for which data, refer to the caption of the figure). We have plotted all 
of the data used, except for F-band, z-band, and K s -band data points, which may overlap 
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with J-band data points if these non-J-band data are converted to J-band. For clarity of 
presentation, we convert non-J-band NIR data to J-band on the basis of the spectral index 
at that point, but still label it with its original band name and plot in a different color. At 
the bottom of the light curve figure, we show, as "residuals," the chi-square contribution of 
each data point. We notice that the radio data provide a large contribution to the total chi- 
square, because of unexpected variations from observation to observation which are difficult 
to reproduce. In particular, the model fails to remain within the several upper limits from 
radio observations (note that eight out of 11 radio observations are upper limits, and three 
are flux measurements). 

Fig. H] also demonstrates the effects of the density on the afterglow evolution (shown 
in dashed and dotted lines, respectively). We notice that the X-ray and optical/NIR can 
be fitted well even for density values varying by 2 orders of magnitude; the only significant 
effects are seen in the radio light curve. The radio lightcurve for a density n = 10 cm -3 
medium peaks around t ~ 2.0 x 10 6 s at a flux of F v m 200 //Jy, that for a density n = 
84.4 cm _3 medium peaks around t ~ 3.0 x 10 6 s at F v lOOyuJy, and that for a density 
n = 10 3 cm -3 medium comes even later, at t ~ 5.0 x 10 6 s and F v ~ 120 /iJy. At this point, 
the afterglow is in the regime with v m < v a < v c , and the peak flux is at the self-absorption 
frequency. Because the self-absorption frequency is a function of density, the larger the 
density, the larger the self-absorption frequency. This predicts an earlier peak for a lower 
density. With more observational data points, or even upper limits after t ~ 3 x 10 6 s, we 
could put stronger constraints on the circumstellar density. 

The posterior distributions for the parameters for model (B) are shown in Fig. [51 As 
for model (A), most of the distributions are satisfyingly Gaussian in shape. One exception is 
the magnetic field equipartition parameter 6B >r , seen to peak at €B, r = 1/2, the upper bound 
for tB, r in our model. Also, we note that the posterior distributions for the density n, the 
initial Lorentz factor r , and the opening half-angle 9, have irregular tails. These irregular 
tail regions correspond to a distinct (local) chi-square minimum. 

In the region where the reduced chi-square for model (B) reaches its minimum value 
53.0/28 « 1.9, the best fit parameters are p = 2.243, e e = 0.0084, e B j = 5.7 x 10~ 3 , 
n = 212 A cm" 3 , 9 = 0.126, E 52 = 146, A v = 0.032 mag, r = 183.6, and e B , r = 0.50. The 
light curve for these parameters is given in Fig. [61 

In Figure [71 we show contour plots for the joint confidence regions of three important 
physical parameters: the jet opening half-angle 9, the isotropic-equivalent kinetic energy 
£752, an d the circumburst density n. This illustrates the degree of covariance between these 
quantities, in a quantitative manner. 
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In Table |3] we list the best-fitting parameters, and the parameter ranges for 1-a and 
90% confidence level for both models (A) and (B). 



4. ANALYSIS OF THE RESULTS 

4.1. J-Band Light-Curve and IC Suppression 

The shift of the afterglow starting time to the epoch t ~ 402 seconds modifies the 
analysis and the interpretation of the early afterglow in model (B), but this should not affect 
the late-time afterglow light curve evolution, which should be similar to the one for model 
(A). Thus the discussion below will be focused on the light curve of model (A), and where 
there are differences, these are pointed out. 



Haislip et al.l (120061 ) have fitted their collection of NIR data, and find that between 
3 hours and 0.5 days after the burst, the fading of the afterglow can be fitted by a power law 
of index a = 1.36^0;^, while after 0.5 days the fading appears to slow down to a temporal 
index of 0.82+^g. At t = 10.6 hours (0.44 days), the spe ctral index is 0„ = 1.2 5+°-^. 



A single power-law decay is ruled out at 3.7a confidence. iTagliaferri et al.l (120051 ) have 
extended datasets whose observation time reaches up to 7 days after the burst, and fit the 
lightcurve with a smoothly-broken power law. The fit gives ol\ = 0.72+Q2Q, ct2 = 2.4 ± 0.4, 
and tb = 2.6 ± 1.0 days. In addition, the spectral index at t = 1.155 days is calculated to 
be (3 — 1.25 ± 0.25 or (3 = 1.2 ± 0.3 by two slightly different fitting codes. Thus, based 
only on the observations in the optical/NIR bands, we can divide the light curve into 3 
segments: (D) t < 0.5 days: the afterglow decays as a power-law with index of 1.36lo;ocs; (E) 
0.5 days < t < 1.6 days: the light curve is relatively flat, decaying as a power-law with index 
a ~ 0.82; and (F) t > 1.6 days: the light curve decays with an index of 2.4 ±0.4 (see Fig. [1]). 



Wei et al.l (120061 ) have argued that the fast decay during stage (D) represents the normal 
afterglow, and the flattening at stage (E) is caused by energy injection. Then the stage (F) 
would indicate a return to the normal afterglow evolution. Here, however, we have presented 
a different interpretation for the stages (D) and (E). The flux at stage (D) is considered to be 
suppressed by the inverse Compton interaction between electrons in the forward shock and 
X-ray flare photons, while the flux in stage (E) is the n ormal flux w i thout external inverse 



Compton process. This is motivated by the argument of I Wang et al.l (120061 ) suggesting that 
the X-ray flare photons can interact with the electrons in the forward shock regions via inverse 
Compton scattering. The origin of the late-time X-ray flares is unknown, although a widely- 
held view is that they are due to the internal shocks from late time engine activity. Since 
these X-ray photons would be coming, in this view, from a region different from (and behind) 
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the forward shock, we call this Inverse Compton (IC) process an "external IC process." In 
this case, the external IC process will contribute significantly to the cooling of the forward 
shock electrons, since the flare luminosity is much larger than the forward shock (afterglow) 
luminosity, and the Compton parameter is determined by the ratio of those two luminosities. 
If the total radiated energy at a given time is constant, when the energy radiated via the 
inverse Compton process increases, the synchrotron radiation should decrease. In effect, this 
can be viewed as the synchrotron radiation having been suppressed by the inverse Compton 
process. Since the J-band luminosity at this time is dominated by synchrotron radiation, the 
observed flux would become smaller in the presence of strong external IC processes. At the 
end of the flare, the external inverse Compton suppression disappears, and the synchrotron 
radiation can then return to its normal course. 

If we assume that the averaged luminosity ratio between the flare and the forward 
shock in the X-ray is k = £ic,fl/Asyn,f ("fl" is the flare and "f" denotes the forward shock), 
following the similar definition for the Compton parameter which is the ratio of the IC lumi- 
nosity to the synchrotron luminosity Issc = ^ic,f/-^syn, f = (MzI§b)]—. wh ere the subscript 



"SSC" indicates the self-Compton scattering process (ISari fc Esinl l200ll ). we can get the 
new Compton parameter, considering the external IC process, as lic,fl = ^ic,totai/-^syn,f — 
[(-^ic.fi + Lic t t )/Asyn,f] = [{k + l)r]e e / es] 1 ^ 2 . We can see an additional factor of (k + l) 1 / 2 
contributes to the Compton parameter for the external IC process compared with the Comp- 
ton parameter for the usual SSC case. Because normally the parameter k 3> 1, we expect 
Yicext > issc- For the fast-cooling case, rj = 1, so we have F ICifl)fast = \(k + l)e e /e B ] 1/2 
and Ysscfast = i e e/ e B) 1 ^ 2 - The IC process will affect the synchrotron radiation, and change 

1 /2 

the cooling frequency v c for synchrotron radiation and the observed flux F v oc v c ■ For 
synchrotron radiation, v c oc (1 + Y)~ 2 and the external IC process will lead to a much lower 
cooling frequency u c due to an increase in the value of the Compton parameter. During the 
time period of stage (D) in Fig. 1, the electrons are in the fast-cooling regime, v c < u m < u a , 
where v Q is the observing frequency, so that the ratio of the flux without external IC to 
the flux with external IC is S = F no /F lc = (1 + Y lCMast )/(l + Y SS cMt) ~ (k + 1) 1/2 , 
where Fjq is the flux with external IC process considered and the flux without IC is F no = 
(l + r ICiCxt )/(l + lssc)^ic,ext ~ (fc + l) 1/2 F Ic . Therefore, the observed flux in the optical/NIR 
band during section (D) should be multiplied by a factor [(1 + k) 1 ^ 2 ] to recover the flux that 
would be observed without IC suppression effects. 

We can estimate the required luminosity ratio from the observed suppression factor. 
Taking the power-law index for the electron energy distribution to be p = 2.2, the theoreti- 
cally expected temporal index is a = (3p — 2)/4 = 1.15 for v > (z/ c , u m ), the expected regime 
for the optical afterglow during stage (E). Then we can extrapolate the optical lightcurve 
back from t = 1.92 x 10 5 s (where the flux is / = 9.14 ± 1.75/iJy) to the observer time 



-20- 



t = 2.7 x 10 4 s, where the observed flux is / ~ 55 /zJy. This theoretical flux, in the absence 
of the external IC process, will be 9.14 ■ (1.92 x 10 5 /2.7 x 10 4 ) 1 15 ~ 81.0 //Jy. Thus, the flux 
will be suppressed by a factor of S ~ 1.6 ~ (1 + k) 1 / 2 . From the observational point of view, 
we invert this problem and solve for k, deriving k « 1.6. 

Apparently, compared to the mean observed luminosity ratio between the flares and the 
forward shock, Ffl/Ff s >10 (see Fig. [TJ, this indicates within our picture that only a small 
fraction of the flare photons have interacted with the electrons in the afterglow. A possible 
explanation could be that the anisotropi c distribution of t he incoming flare photons in the 



comoving frame of the afterglow shock (jWang et al.ll2006l ). This resultsing more head-on 



scattering which reduces the IC interaction. Therefore, the suppression is relatively small 
and the optical flux is only affected to a reduce degree. 



4.2. Radio Light-curve 



GRB 050904 shows several similarities with GRB 990123, including a large isotropic- 
equivalent gamma-ray energy and a very bright optical flash. In GRB 990123, the radio 
emission was observed to peak at t — 1 day in the observ er frame, and this emission was 
interpreted as the radio emission from the reverse shock (ISari &: Piranl Il999l ). Given the 
important role of the reverse shock in our model (B) of GRB 050904, it is interesting to 
consider its radio emission. However, our estimates indicate that the reverse shock radio 
emission at early time would be suppressed by the large circumburst density in our model, 
so that we would not expect to observe a bright radio flare. In Fig. [6l we have plotted the 
reverse shock radio emission as a dotted blue line starting from t = 470 s, while the solid 
blue line shows the combined radio emission from both the forward shock and the reverse 
shock; as can be seen, the emission from the reverse shock is negligible. 

In Fig. [61 there is one averaged radio data point centered at t ~ 2 x 10 5 seconds and 
/„ ~ 30 [A Jy (indicated with dash ed lines ) based on the average of multiple VLA observations 
spanning ~ 6 days (see Fig. 2 of iFrail et al.ll2006l ). Given the large dynamic range in time, 
we consider this data point as an upper limit of the radio emission during that period. We 
can see that our best-fit radio light curve for the forward shock in both models (A) and (B) 
fits accommodates this upper limit reasonably well, along with the other data points. 
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4.3. Density and Energy Constraints 



The density obtained in our fit, n ss 84.4 cm 3 for model (A) an d n ~ 212 . 4 cm 3 
for model (B), is smaller than the n = 680 cm -3 density derived by iFrail et al. r|20Q6h . 
Investigating our model fits, we find that the discrepancy arises because their model peaks 
later in the radio than ours. Referring to the posterior distribution for the density parameter 
in Figs. [3] and [5J we find that the most likely range (90 %-confidence) for density is from 26 
to 273.4 cm -3 for model (A) and from 87.8 to 270.6 cm -3 for model (B). On the other hand, 
the light curves for different densities demonstrate that there is no big noticeable affect of 
these density changes on the X-ray and optical/NIR lightcurves. 

The radio observations thus provide the best frequency for density constraints, and fail 
to provide a tight constraint mainly because of the lack of data (measurements) late times. 

With regards to the circumburst density, we note that iKawai et al. observed 
the spectrum and detected in it fine structure lines including Sill*, which were taken to 
imply an electron density of up to io 2 - 3±0 - 7 cm~ 3 . The density obtained from the spectral 
lines would thus be consistent with our best-fit value for the density, and suggest that the 
observed lines formed in re gion similar to tha t hosting the GRB itse lf. However, we note 
that several recent papers ( Berger et al. 2005 ; Prochaska et ali bood ) have suggested that 
these fine-structure transitions in GRB afterglows are excited by radiative processes, rather 
than collisions, which would make the density constraint irrelevant. 



An important check on our model results, pointed out by IFrail et al.l (120061 ). is to 
estimate the X-ray luminosity from the X-ray light curve at some fiducial time (usually 
at t — 10 hours). At t — 3.36 days, the observed flux over the XRT energy range is 
~ 2.1 x 10~ 14 erg s _1 cm -2 . Calculating back to t ~ 10 hours for the afterglow flux from 
synchrotron radiation, the predicted flux is / ~ 2.8 x 10~ 13 erg s" 1 cm" 2 . Therefore the 
X-ray luminosity should be L x ,i SO = 47r£)|/(l + ^)- Q +/ 3 - 1 ~ 1.66 x 10 46 erg s -1 , and the 
geometrically-corrected X-ray luminosity we find is different from and significantly lower 
than theirs. The reason for this is that when they calculate the isotropic-e quivalent en 



ergy, they appear not to have included the /c-correction factor (1 + z)~ a+l3 ~ l . iBerger et al. 
(120031 ) also appear not to have included the ^-correction factor when doing statistics on the 
isotropic-equivalent and geometrically-corrected X-ray luminosity. Here we do the statistics 
again after putting the correction, and find that the geometrically-corrected X-ray luminosi- 
ties are still clustered, but the peak value has shifted down to L x , p ~ 10 44 erg/s, reduced by 
a factor of five. 



Once we have obtained the X-ray luminosity, we can estimate the kinetic energy with 
the best fitting parameters and the observed spectral index. In model (A), we have best 
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fitting parameters e e = 0.0309, e B j = 0.198 , /3 ~ 1, p = 2.152, and a = (3p - 2)/4. From 
the equation ID4I in the appendix, we have E 52 — 24.4, which is consistent with the kinetic 
energy derived from the fit, Ek = 22.4. This is not surprising since the formula in the 
appendix is only a shortcut to obtain the kinetic energy. In model (B), we have the best 
fitting parameters e e = 8.4 x 10 -3 , e B j = 5.7 x 10 -3 , /3 ~ I, p = 2.243, and a = (3p-2)/4, 
and similarly we obtain Em — 148. After considering the fc-correction factor for the X-ray 
luminosity by iFrail et al.l (120061 ). we can re-estimate the kinetic energy as Em — 230 with 
the parameters e e = 0.02, and e B = 0.015, which is 3 times larger than their best-fit kinetic 
energy E 52 — 88. 

Considering our estimated opening half-angles of 9 ~ 0.13 for both models (A) and 
(B), we can calculate the geometrically-corrected X-ray luminosity for GRB 050904 to be 
Lx = Lx,iso x 2 /2 « 1.40 x 10 44 erg s~\ which falls within the corrected luminosity range 
of low-redshift GRBs faerger et alJl2003h . 



We now calculate the geometrically-corrected kinetic energy. Our best-fit opening half- 
angle is 6 ~ 0.13 for both models, and the fitted kinetic energy for model (A) and (B) is 
2.24 x 10 53 ergs and 1.47 x 10 54 ergs, so the geometrically-corrected kinetic energy will be 
1.84 x 10 51 ergs and 1.17x 10 52 ergs for models (A) and (B), respectively. Broadband modeling 
of 10 low-redshift bursts indicated that the geometrically-corrected kinetic energies of two 
were anomal ously high, 2 x 10 51 ergs, a pproximately 10 times higher than for the other 
eight GRBs (IPanaitescu fc Kumarl |2002| ) . Our geometrically-corrected kinetic energy for 
model (A) is comparable to those anomalously-large kinetic energies seen from low-redshift 
GRBs. The kinetic energy of model (B), on the other hand, is significantly larger even than 
this. Both models yield a relatively large kinetic energy has been obtained for GRB 050904. 
GRB 050904 thus suggests that bursts at high redshift are somehow able to tap into a higher- 
energy reservoir than the low-redshift events. 



4.4. Burst Energetics and Efficiency 



The radiated isotro pic-equivalent gamma-r ay energy for GRB 050904 is 6.6 x 10 



53 



< 



£ 7iiso < 3.2 x 10 54 



ergs (ICusumano et al.ll2006l ). If the isotropic-equivalent kinetic energy 



for the afterglow of GRB 050904, as we concluded in case (A), is 2.24 x 10 53 ergs, we ca n 



estimate the GRB efficiency as ( = E^^ SOj52 / (E 1)isQj52 + E 52 ) (ILloyd-Ronning fc Zhang!l2004j ) 



which is 74.7% < ( < 93.5% for GRB 050904. If the isotropic-equivalent kinetic energy of the 
afterglow, as in case (B), is 1.47 x 10 54 ergs, then the corresponding GRB efficiency is roughly 
31.0% < C < 68.5%. In either case, this indic ates that GRB 050904 ha s a hig h efficiency; 
however, such high efficiencies are not unique. ILloyd-Ronning Zhang! (120041 ) found that 
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a substantial number of GRBs have high efficiency. Fo r some bursts like GRB 990705, the 
inferred efficiency even reaches up to 99%. iFan &: Piranl (120061 ) have shown that the inferred 
efficiency can be r educed when inverse Compton effects are taken into account (see also 



Granot et al.l 120061 ) . Even so, there are still several bursts which ha ve a high observe d 



efficiency, for example, 94% for GRB 050315 and 80% for GRB 050416 (IZhang et all 120070 . 



5. Discussion 



The extreme interest in GRB 050904 has motivated several groups to analyze the burst 
data and suggest interpretations. These works fall into two categories: (1) Comparing the 
properties of GRB 050904 with other bursts; and (2) Making fits to the GRB 050904 data 
either in the framework of internal shock or external shock models. We give a brief description 
below of the work of other groups and mention key differences between their work and ours. 



In category (1). iKann et al.l (120071 ) made a composite J-band light curve starting from 
2 x 10~ 3 days to ~ 23 days. After applying extinction correction, they shifted GRB 050904 as 
well as lower-redshift bursts to z — 1, and made a comparison of afterglow lightcurves. They 
found that GRB 050904 is much brighter than other GRBs at early times, but of roughly 
equal brightness at late times. Thus they conclude that GRB 050904 most likely is still a 
normal GRB. 



The other analyses are all in category (2). IZou et al.l (120061 ) argued that GRB 050904 is 
a burst with extremely long central engine activity. They put all the observed data within 
the framework of the internal shock model. By contrast, we only treat the first several 
hundred seconds (BAT) as internal shock activity in our model (corresponding to the stages 
(A) and (B) for model (A), and stage (A) for model (B) in Fig. 1). The late-time X-ray flares 
between 6 x 10 3 and 6 x 10 4 seconds may be due to internal shock, but we have interpreted 
portions of this X-ray emission as being due to the forward-shock afterglow. 



Wei et all (120061 ) argued that the t 



470 s flare is from internal shocks on the basis of 
its fast decay (the temporal index is a ~ 8.8 relative to the trigger time) and also because 
the opt ical-to-X-ray em ission of the flare cannot be described by a synchrotron radiation 
model (jBoer et al.ll2006l ). They made fits to all the available J-band data. They argue that 
the slow-decay portion of the lightcurve is due to energy injection. In our model, besides 
fitting over all the additional bands (X-ray and radio), we propose a new mechanism for 
the flattening, namely that it is caused by the suppression of the synchrotron radiation by 
the interaction between the X-ray flare photons and afterglow electrons. Separately, in our 
model (B), we introduced a new reference time t' which flattens the decay index, and allows 
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us to interpret the t fa 470 s optical/X-ray flare as arising in the reverse shock. 



Frail et al.l (120061 ) made broadband model fits including the X-ray, NIR/optical, and 
radio data. The difference with our work is partly that we have used the larger data set that 
later became available. We have included two X-ray data points as early as t ~ 10 3 seconds. 



and w e also included the H band data observed t ~ 23 days after the burst (IBerger et al. 
20061 ). Since most of the other data are concentrated around t ~ 10 5 seconds, the introduc- 



tion of these additional X-ray and NIR data have some impact on the final fitting result. 
The other major difference is that we have freed all the possible parameters and applied the 
MCMC method for the global fitting, making an efficient exploration of the full parame- 
ter space, and providing the posterior distributions (including confidence intervals) for each 
parameter. 



Gendre et al.l (120071 ) argued that the power-law-like decay right after the flare t ~ 470 s 
should be interpreted as forward shock emission. Because the extrapolation of the late time 
X-ray flux to early times is lower than the observed value, they found that a wind-type 
environment was favored by the closure relationshi p for this early-time seg ment (it should 



also be noted that their spectral index is smaller than lCusumano et al.l (120061 )). They propose 
a density-jump model for the afterglow evolution: before a certain radius, the density goes 
as n oc r~ 2 and after that, the density is a constant ISM model. At the transition point a 
termination shock is formed which lies around Rt ~ 1.8 x 10~ 2 pc from the central engine. In 
our model, we consider the same segment of data, but interpret it differently. We argue that 
the flux between 600 and 800 seconds arises from the combination of three sources: high- 
latitude emission of the prompt emission, flux from the flare, and forward shock emission. 
Reviewing the data closely, we see that actually there are two other small flares between 
800 and 2000 seconds, so the data around t ~ 1500 seconds has the contribution from the 
flares and the forward shock (see Fig. 2). Once we subtract the flare contribution from 
the observed data, we argue, the forward shock contribution is what remains. And in fact, 
we find that an extrapolation of the late-time flux to early times is consistent with this 
flare-subtracted early-time flux. 

There remain some substantial differences in parameter values between the two models 
that we have presented. For example, the most likely value for esj in model (A) is roughly 
~ 0.1, but the most likely value for ebj in model (B) is around 5 x 10~ 3 . It turns out that 
the main difference lies in the optical flux at early times. We take as an example the J-band 
light curve for the best fitting parameter set in model (A). The flux slowly decays with an 
index of a « 0.5 before t ~ 2 x 10 3 and then follows a faster decay with a temporal index 
of a ~ 1.1 after t ~ 2 x 10 3 . The break at t ~ 2 x 10 3 seconds is caused by the crossing of 
the electron's typical frequency v m through the optical observing frequency. Since we have 
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1/2 

v m oc e B j, then the smaller sbj, the smaller the typical frequency, therefore the earlier the 
crossing time. The fit in model (B) requires an earlier break than in model (A), so a smaller 
value of tBj is expected in model (B). Since the observed flux from the forward shock is 
the same in both models, we expect a higher kinetic energy in model (B). Similarly, at the 
deceleration time, the reverse shock flux is larger than that from the forward shock, so a 
larger ee, r is expected in the reverse shock. 

We notice that in both our models we find a small value for e e , e.g., the best fitting 
e e = 3.1 x 1CT 2 for model (A) and e e = 8.4 x 1CT 3 for model (B). Considering the radiative 
correction factor R = (t/t ) 17ee ^ i8 (where t is the deceleration time, and we set t = 300 
seconds for model (A) and to — 470 seconds for model (B)), that is, the amount by which 
the kinetic energy at time t is reduced by comparison to the kinetic energy at time to, we 
find R = 1.07 for model (A) and R = 1.02 for model (B) at t = 1.0 x 10 5 seconds after the 
burst. Therefore, the radiative losses are mild in either case, consistent with our assumption 
that the afterglow evolves in an adiabatic fashion. 

Our fitting results also show that the host galaxy dus t extinction is quite small, Ay ~ 0.1 



mag or even smaller, consistent with the results from iKann et al.l (120071 ). They applied 
several different dust models (MW, LMC and SMC) over the composite J-band, and all 
the models suggested zero or negligible extinction. In the context of GRB 990123, it was 
suggested that the negligible extinction to th e burst was the re sult of dust destruction by 



the strong burst and early-afterglow emission (IKann et al.ll2006l ). 



6. CONCLUSIONS 

In this paper we have performed the most extensive multiband analysis so far of the 
GRB 050904 afterglow. We have considered two scenarios: (A) Only forward shock emission 
is considered, and the flares peaking at t ~ 470 after the burst are assumed to be due to 
internal shocks (or are otherwise independent of the afterglow); and (B) The NIR and X-ray 
flares at t « 470 s are ascribed to emission from the reverse shock - when the ejecta has 
swept up enough material and starts to decelerate, the synchrotron radiation in the reverse 
shock produces the optical flare, and the self-Compton scattering of synchrotron photons 
generates the flare observed in the XRT and BAT energy bands. 

Combining the early afterglow data with late-time observations in the X-ray, optical 
and radio, and using a Markov chain Monte Carlo method, we present a full characterization 
of the posterior distributions (including confidence intervals) for the various parameters of 
our model fits. Our best-fit parameter values for model (A) are p = 2.152, e e = 0.309, 
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e BJ = 0.198, n = 84.4 cm" 3 , 9 = 0.128, E 52 = 22.4, and Ay = 0.0343 mag, with a 
reduced cm-squared value of 36.2/26 ~ 1.39. Our best-fit parameter values for model (B) 
are p = 2.243, e e = 0.0084, e BJ = 5.7 x 10~ 3 , n = 212.4 cm" 3 , 9 = 0.126, E 52 = 147, 
Ay = 3.18 x 10~ 2 mag, T = 183.6, and e B , r = 0.50, with a reduced chi-squared value of 
53.0/28 ~ 1.89. Note that the subscripts r and / refer to the reverse shock and forward 
shock respectively, e B j = e B}r /R B , anc ^ we have assumed e e j = e e>r . 

We have compared the density, the geometrically-corrected kinetic energy and the X-ray 
luminosity at t — 10 hours derived here for our two models of GRB 050904 against those 
values for other bursts, as derived from afterglow modeling. The results for both models 
show that although the X-ray luminosity of GRB 050904 falls within the range for low- 
redshift GRBs, the density and geometrically-corrected kinetic energy are both above the 
typical values for low-redshift GRBs, which suggests that GRB 050904 may be a member of 
a distinct population of high-redshift, higher kinetic-energy bursts, whose properties differ 
from those of low-redshift GRBs. A clear preference between our (A) and (B) models is 
hard to establish at present, since there is only one high-redshift GRB known. One would 
like access to several more high-redshift GRB observational datasets before attempting to 
discriminate between the two models. 



It is estimated that ~ 7 % — 40% GRBs are located at z > 5 (IJakobsson et al.ll2006l ). and 
detection rate simulations bv lGou et al.l (120041) indicate th at Swift could detect GRBs out to 
redshift z ~ 30, if they are present. iBromm &: Loebl (120061 ) also predict that 10% of the Swift 
GRBs originate at z > 5. It appears that one can realistically expect a handful (5 to 10) 
of additional high-redshift GRB detections with rapid follow-up in the next few years of the 
Swift mission. In this case, the consistent application of MCMC methods, as used here, will 
lead efficiently to a set of statistically well-quantified, posterior parameter distributions and 
confidence intervals. This would enable a statistically meaningful comparison of high-redshift 
and low-redshift GRB parameters, which might well lead us to a definite understanding of 
the physics and environments of GRBs as a function of redshift, up to the highest redshifts 
detected. This would also have a substantial impact on the study of the large scale structure 
and star formation processes throughout the Universe, and the properties of the cosmic 
reionization at z ~ 6. 
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A. Self- Absorption Frequency 

Once we consider radio emission in our afterglow models, the self- absorption frequency 
becomes an important parameter to consider. Below we give the expressions for the self- 
absorption frequency in the different regimes. 

If we assume that the electron distribution follows a power law A/"( 7e ) = A^ xT p where 



7i < 7 e < 72, then the self-absorption coefficient for the various possible regimes is ( jWu et al. 
2003h : 



-(P+4) 



W p+4) (-)~ (P+4)/2 (Al) 

-(P+4) ( u_Y h ' 2 p -v/»2 



-5/3 



k = |iv 7 < 



where Cl = ^^^Sts* <* = ^^(p + f W)^), c 3 = ^(p + 2), v x 
and ^2 are the typical synchrotron frequencies of electrons with the Lorentz factor ji and 
72, respectively, and T(x) is the Gamma function. 

Following the definition of K v {v a )l = 1 for the self-absorption frequency (where I = R/Y 
is the thickness of the shell), we can find the self- absorption for the forward shock region. 
The self-absorption frequency for the reverse shock has a similar form, with the difference 
that the quantities specific to the forward shock region should be replaced with those specific 
to the reverse shock. 

Fast Cooling: v a < u c < u m , 

•P-oJ^ggT (A2) 



3B 7c 5 / 

where the superscript "i" denotes the different regimes 
Fast Cooling: v c < u a < u m , 



= ^(^) 1/3 (A3) 



c 2 q e noR 

^ = ^ /9 W ) ) B/9 (A4) 



Fast Cooling: v c < u m < u, 



a ■ 



/ in 



& = Vm | =p , (A5) 



This preprint was prepared with the AAS IATgX macros v5.2. 
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^ = [(^ 1) ) 10/3 ^ 8/3 <- 1 ] 1/(p+5) (A6) 
Slow Cooling: v a < v m < u c 

We notice that v\lv\ — (p — l) 3 / 5 (i/ c /z/ m ) 1 / 2 . To keep the continuity of the flux while the 
afterglow transits from the fast-cooling regime to the slow-cooling regime, we divide vt^ by 
a factor of (p — 1) 3//5 . 

Slow Cooling: v m < v a < u c , 

2/(p+4) 



^Ht) (A8) 
«f = My%^v<^>» <A9) 



Slow Cooling: v m < v c < u a , 

,«,_„ f c,(p-l) qe n R Ym ~^ 2/{P+5) (A1Q) 



3£? 7 r 4 

V f) = [(z.( 1 )) 1 °/3 zy 8/3 zy p-l ] l/(p + 5) (An) 

it can be shown that the self-absorption frequency in the regime v m < u c < v a has the same 
form as that for the fast-cooling case, u c < v m < u a . 

Because the observer time, t obs , at z = is connected to the time in the source frame, 
t s , at redshift z by the relation t b s = (1 + z)t s , we have the redshift dependence for main 
characteristic quantities: the shock radius R oc t^ 4 oc (1 + <z)~ 1//4 , the shock Lorentz factor 
T oc tj 3 ^ 8 oc (1 + z) 3 / 8 , the magnetic field B oc T oc (1 + z) 3//8 , the typical Lorentz factor 
7 m oc T oc (1 + z) 3 / 8 , the cooling Lorentz factor 7 C oc t^T -3 oc (1 + z)~ l l 8 . In addition, we 
have v m oc (1 + z) 1 / 2 and v c oc (1 + z)~ x l 2 (see Eqn. [2]). Substituting these dependence into 
the relations for the self-absorptions above, the redshift dependence for the self-absorptions 

is OC (1 + Z)- 1 / 2 ,^ OC (1 + Z)- 1 / 2 ^ 3 OC (1 + z )b-7)/[2(p+5)] )Z ,(4) ^ (1 + ^-1^(5) ^ 

(1 + z )(p-6)/[2( P+ 4)] > and u f) = u f) K (! + 2 )b-7)/[2(p+5)]_ 



B. Inverse Compton Spectrum 



As described by ISari &: Esinl (1200 ll ) , the inverse Compton flux can be calculated from 
the following double integral: 



f 



IO 



Ext. 



T 



djN^) / dxf v ,{x) 
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where N(j) is the electron distribution in the shocked shell, f Va { x ) is the seed photon flux, 
and x is defined as x = v/A^ 2 v s where the subscript "s" denotes the seed photon. 



B.l. v c < u a < u m 

The distribution of seed photons is described by the synchrotron spectrum, a broken 
power law with the characteristic quantities (Sari et al. 1998). Then the inner integral in 
Eqn. fIBll) gives: 



A — ■j/mm^O ( t o (21 



1 2 — ^fmax-EO 
2 



4-y 2 u m xo 

^ 3 ~ (p+2) fmax%0 Cffi J (^472 



47 2 z/^ 2) a;o < 1/ < 47 2 z/ m x 
z/ > 47 2 i/ m x . 



:B2) 



The integration over different electron energies again needs to be divided into four different 
regimes: 



x < 



djN(-f)h+, 
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Evaluating the integrals in Eqn. (1B3I) . we only keep the dominant terms: 
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B.2. v c < v m < v a 
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It is noted that a factor of l/(p — 1) has been multiplied for the slow-cooling case in 
order to keep the IC flux continuous while the afterglow changes from the fast-cooling regime 
to the slow-cooling regime. 



C. Derivation of Radiative Correction Factor 



In the observer frame the ene rgy loss rate is equa l to the rate at which energy is supplied 
to the unshocked matter 4nR 2 p2 (ICohen et al.lll998l ). multiplied by e, the fraction of energy 
that each particle has lost (assuming that the relation between the radius and the Lorentz 
factor is described by R = AT 2 ct where A=2 for the GRB prompt phase and A e [3, 7] for 
the afterglow deceleration phase). 



dE 



-AnR l p 2 t = -4vrA 2 r 4 c 2 t 2 p 2 e 



(CI) 



where R is the radius of the shocked region from the center in the observer frame, p 2 = 
cT 2 U' /3 is the pressure in the shocked region in the observer frame, and U' = 2T 2 nm p c 2 is 
the energy density of the shocked region in the comoving frame. 

Normally we assume that all the energy stored in elec t rons has been radiated during 
the fast cooling phase, i.e., e = e e . However, I Cohen et al.l (119981 ) have shown that only a 
portion of the energy will be lost even if the electrons are in the fast-cooling regime (refer 
to their Fig. 6). The relation between the radiation factor e and the electron equipartition 
parameter e e is given by e e = 1 — ( ^j= ) 7 (their Eq. 46) where 7 = 4/3 for the extreme 
relativistic case is the adiabatic index. We can solve the above equation numerically, e.g., 
e = 0.11 for e e = 0.2, e = 0.0255 for e e = 0.05, and e = 0.005 for e e = 0.01. We find that 
for the reasonable range of parameter values, we have the following relation, e ~ |e e , and 
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substituting it and BM solution back into Eq. IC1I above, we obtain: 

dE 17 E 

~dt ~l2A te J 



(C2) 



We have the solution as E{t) = E (-^)~ 17ee ^ 12A where to is the deceleration time. Be- 
cause we have applied the Blandford-McKee solution during the whole derivation, the above 
solution is only valid after the deceleration time. Normally during the afterglow deceleration 
phase, the we set the parameter A = 4, and thus we obtain the radiative correction factor 
after prompt phase R = (t/to) 17{Ee//48 - We note that the above correction factor can only be 
applied to the fast cooling case. Taking e e = 0.1 and to = 300 s, the correction factor at 
observer time t = 10 hours would be R ~ 1.18. It should be mentioned that if the estimated 
correction factor is much larger than the unity, one should not estimate the kinetic energy 
only by introducing the radiative correction factor R for the adiabatic solution, and should 
reconsider the whole evolution including the radiation loss effects. Because at this point the 
radiation is not negligi ble any more, th e afterglow evolu tion will deviate substantially from 
the BM solution (e.g., hfost et alJliooSlWu et aliboosh . 



The radiative correction factor R = (t/tn ) 17ee ' 48 for the afterglow is different from the 
one, R = (t/t ) 17ee/16 , provided by ISaril (119971 ) (see also iLlovd-Ronning fc Zhang! I2Q04J ) for 
two reasons: (1) They have taken the value A = 16, while we have used A = 4 which is 
believed to be mu ch more reason a ble from the detailed studies o f the hydrodynamic evolution 
of the afterglow (jWaxmanl 1 19971 : iPanaitescu &: Meszaroslll998f ). (2) They have assumed all 
the electron energ y will be radiated a way if the afterglow is in the fast-cooling regime, but 
from the results of I Cohen et al.l (119981 ). they showed that roughly half of the electron energy 
is lost. 



D. Kinetic Energy 

If the X-ray band is above the typical frequency u m and the cooling frequency u c , the 
observed flux at a certain observer time is given by: 

F = p (p-l)/2 1/2 -p/2 

= 2( 2 " 3p )/ 4 x 10- 30 ergs-W^z" 1 

x(l + ^^^^^r^'e^i^ 74 ^^^ 74 ^^! + Y)-* Dl % (Dl) 

where p is the energy distribution index, E52 is the isotropic-equivalent kinetic energy, e e is 
the electron equipartition parameter, €3 is the magnetic field equipartition parameter, t is 
the observer time, v is the X-ray observing frequency, and Y is the Compton parameter. 
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The X-ray luminosity is described as (ILamb fc Reichartl I200CH : iBerger et al.l 120031 ) : 



L v>x = 4nD 2 L F VtX {l + z )-+^ (D2) 

and 

L x = J L v>x dv 
= CuL UjX 

« 2.5 x 10 45 erg s" 1 (1 + z )(p/^-+P) 

(p+2)/4 (p-1) (p-2)/4 ,(2-3p)/4 (2-p)/2/, , V \-l / m \ 
X - & 52 e e -1 e B,-2 r iOh ^18 U + r J l^J 

where C = J L v ^ x dv jvL v ^ x is an integration constant, and it shows that C ~ 4.3. 

We can find the kinetic energy E from the above equation reversely, and considering 
the radiative correction factor, we obtain, 

E K = 10 52 ergsi?[ ^ f(P+Vh + z) 4(p/4-l/2-a+«/( p+ 2) 

2.5 x 10 45 ergs s _i 

xe 4(p-l)/(p+2) e -(p-2)/(p+2) t (2-3p)/(p+2)^2-p)/(p+2)^ 1 + y^4/(p+2) 

where i? = [^]( 17 / 48 ) £e is the radiative correction factor, and the derivation is given in 
Appendix O 
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Fig. 1.— The combined light curves of GRB 050904 in the BAT, XRT, J and I bands, in the 
observer frame. The B AT (violet, empty triangle) and XRT (red, empty circle) data are taken from 



Cusumano et al.1 ( 2000 ). the early time /-band flare at 9500 A (orange, sta r) is from 



Boer et al 



(2006), and the J-band data (green, solid circles) are from lHaislip et al.l (|2006l ) and lTagliaferri et al 



20051 ). For ease of presentation, we show the BAT flux density divided by a factor 100. The black 
arrows at t = 398 seconds after the burst indicate the adopted reference time point, in our model 
(B), for the start of the afterglow evolution. The solid line (green) shows the observed afterglow 
evolution in J-band. The dashed lines show the theoretically-expected afterglow evolution in the 
X-ray and J-band, respectively, and the dash line (green) is the light curve behavior in the J band 
without the external IC process by the X-ray photons from the flare. The light curve is divided into 
6 sections labeled from A to F sequentially: (A) t < 225 s: prompt emission; (B)225 < t < 600 s: 
simultaneous flares in the X-ray and NIR bands; (C) 600 < t < 6220 s: power-law decay in the 
X-ray band; (D) 6220 < t < 4.3 x 10 4 s: energetic X-ray flaring activity; (E) 0.5days < t < 1.6days: 
flattening of the light curve in the NIR band; (F) t > 1.6 days: a jet break is apparent. Approximate 
temporal power-law indices (F u ~ t~ a ) in the optical (a ) are noted. Estimated spectral indices 
(F v ~ v~P) in the optical/NIR (observed and uncorrected for extinction) and X-ray are also 
provided, whe re available, in the l ower figure panels below the figure. References for the spec tral 
indices are (a) Irlaislip et all j^QQfl); (b) iTagliaferri et al.1 iboOfih: fcUfllCusuma.no et al.1 ifcoOfih . A 
jet break time of ij e t = 2.6 ± 1.0 days has been reported by ITagliaferri et al.l ()2005l ): however our fit 
of all available data suggests instead tj e t = 3.17 ± 0.22 days. 
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Table 1: Reverse and Forward Shock Observation Data Points 



No. 


Obs. Time 


Obs. Freq. 


Band 


Flux(Error) 




logl0(t[sec.]) 


loglOfz/fHzl) 

v o v \ L J/ 




( uJv) 


it 


2.66 


14.6 


I 


5 idol xlO 4 




2.66 


18.1 


5 keV 


70.4 (44) 


:; 


2.66 


19.1 


50 keV 


4 45 (89) 

J- • JI KJ 1 Ul/ / 


4t 


2.86 


18.1 


5 keV 


65(26) 

KJ • KJ KJ \ £J KJ / 


5 f 


3.25 


18.1 




0.32(7) 


6 f 


4 75 


18 1 

-L KJ . ± 




32(40) x10~ 2 

\J .KJ^j \ ±\J I A 1U 


V 


5.46 


18.1 




0.51(17)xl0~ 3 


8 


2.78 


14.6 


I 


12.8(64) xlO 3 1 


9* 


4.98 


14.4 


J 


19.1(73) 


10 


5.27 


14.4 




9.2(17) 


11 


5.56 


14.4 




2.74(19) 


12 


5.66 


14.4 




1 67C32) 

XiW 1 \ KJ J 


13 


5.79 


14.4 




0.42(21) 1 


14 


4.98 


14.27 


H 


23 70 (21) 

KJ ■ 1 KJ \ _L / 


15 


5.26 


14.27 




10 60 (40) 


16 


6.30 


14.27 




8 0(25) xlO -2 

KJ • \J \ _ KJ J s\ J. \J 


17 


4.99 


14.14 


± -5 


33 70(210) 

KJKJ • 1 U 1 1 


18 


5.27 


14.14 




15 0(10) 


19 


4.94 


14.51 




9 12(19) 


20 


5.03 


14.51 




6.92 (10) 


21 


4.97 


14.46 


Y 


14.0(30) 


22 


5.29 


14.46 




8.4(22) 


23* 


4.64 


9.9 


Radio 


89.0(580) 


24 


5.08 


9.9 




41.0(250) 


25 


5.67 


9.9 




-3.0(250) 


26 


5.73 


9.9 




27.0(240) 


27 


6.24 


9.9 




89 0(370) 

KJ t_J • \J \ KJ 1 KJ J 


28 


6.40 


9.9 




40.0(300) 


29 


6.46 


9.9 




-10.0(350) 


30 


6.47 


9.9 




64.0(230) 


31 


6.48 


9.9 




116.0(180) 


32 


6.51 


9.9 




67.0(170) 


33 


6.58 


9.9 




13.0(270) 



iThe reverse sho c k emi ssion. The NIR data are from Boer et al. ( 200(jh . and the X-ray data are from 
Cusumano et al. ( 20061 ). 



^The X-ray data are from Cusumano et al. (2006). For the early-time X-ray data, the contribution from the 
flares has been subtracted. 

*The J-band data are from lHaislip et aL ( 2006 ') and Tagliaferri et al. ( 20051 ). Groups of adjacent data points 
have been averaged. 



*Radio data are from lFrail et al.l ( 2006T ). 

^Assuming both upper limits as 2-a limits, we convert the upper limits into synthetic measurements with 
error bars by taking the measurement to be half the upper limit, and the error bar to be one-quarter the 
unnfir limit. 
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Fig. 2. — Illustration of the multiple mechanisms contributing to the observed flux at early times. 
The high-latitude component of the prompt emission is described by a broken power-law with a 
break time of t « 350 s in the observer frame (jCusumano et al.ll2006l ). Before the break time, a ~ 2 
and afterwards, a ~ 3. The main flare peaks at t = 468.0 ± 2.0 s after the burst. The black data 
points correspond to the observed flux. The green lines indicate the theoretical extrapolation of the 
high-latitude emission to late times, using the best-fit values of the burst and flare emission without 
any subsequent fitting. The red-color boxes indicate the data used for estimating the forward shock 



emission. 
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Table 2: Other observational constraints for the fitting 





Constraint 


Comments 




Reference 




= 3.17 ±0.22 


Jet break time 




1 


^dec 


= 468.0±2.0 2 


Deceleration time 




2 




= 1.25±0.25 3 


Spectral index in J band at t= 


1.155 days 


3 




= 0.96±0.19 4 


Spectral index in X-ray from t 


=680 and 1600 s. 


4 



References. — ( 1) The jet break is derived from fitting over all available data in the X-ray and NIR/optical. 
Kann et all |2007h find t jet = 2. 63±0.37 days (1-a), wh ile ITagliaferri et all (Eoqih fou nd < jct = 2. 6 ±1.0 days. 
(2) See the text in Sec. EH (3) ITagliaferri et all |2005l ). (4) ICusumano et all l|2006| ). 



Table 3: The best fit values and parameter ranges for models (A) and (B). 



Parameters 


(A) Forward shock only 


(B) Reverse shock flare 


Best Fit 


1 a Range 


90% Range 


Best Fit 


1 a Range 


90% Range 


V 


2.15 


2.11-2.19 


2.09-2.22 


2.24 


2.20-2.29 


2.18-2.32 


ee (/io- 2 ) 


3.09 


4.3-14.6 


2.8-26.3 


0.84 


0.75-1.3 


0.66-1.6 


esA/lO- 2 ) 


19.8 


4.5-38.9 


2.0-50.5 


0.57 


0.32-0.58 


0.26-0.70 


n 


84.4 


26-273 


9 - 580 


212.4 


88-271 


58-470 


e 


0.128 


0.12-0.18 


0.10-0.19 


0.126 


0.11-0.13 


0.11-0.14 


E52 


22.4 


13-53 


7-102 


146.6 


114-182 


93-208 


Av(/10- 2 ) 


3.43 


1.8-8.0 


0.7-10.6 


3.18 


1.9-7.8 


0.7-9.6 


r 








183.6 


176-206 


163-219 










0.50 


0.4-0.5 


0.3-0.5 


X 2 /dof 


36.2/26 






53.0/28 
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Fig. 3. — Posterior distribution of all parameters for model (A), which excludes data related to 
the flaring activity at t ~ 470 s: p, the power-law index of the electron energy distribution; e e , the 
electron equipartition parameter; the magnetic field equipartition parameter in the forward 

shock; n, the circumstellar density in units of cm -3 ; 6, the opening half-angle before jet break; 
#52, the isotropic-equivalent kinetic energy in units of 10 52 ergs; and Ay, the dust extinction of 
the host galaxy. In each plot, the shaded blue region delimits the \-a (68.2%) confidence range, 
and the vertical lines indicate the 90%-confidence range. The green color indicates the posterior 
distribution of the parameters for models having j m < 2.1 at t = 10 7 s, with the red lines indicating 
1-cr confidence ranges, and the height of the distribution magnified by a factor of five for clarity; 
these model realizations correspond to roughly 1.5% of the total. 
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Fig. 4. — Top: Theoretical light curves (solid lines) corresponding to the best-fit parameters for 
model (A), which excludes data related to the flaring activity at t ~ 470 s. Best-fit model parameters 
are j> = 2.152, e e = 0.031, e B ,f = 0.198, n = 84.4cm" 3 , 9 = 0.128, E 52 = 22.4, Ay = 0.0343 mag. We 
have shown all the available data on the plot. If the data point is used in the fitting, it is plotted 
with a bright color, otherwise it is plotted in grey. To illustrate the effects of different densities for 
model (A), we show additional light curves. The dashed lines correspond to n = 10 3 cm -3 , with 
marginalized best-fit parameters of p = 2.25, e e = 0.019, e^j = 0.043,0 = 0.182, E 52 = 44, Ay = 
0.0426 mag, and with reduced chi-squared value 46.1/26 = 1.77. The dotted lines correspond to 
n = 10 cm -3 , with marginalized best-fit parameters of p = 2.19, e e = 0.014, e#,/ = 0.078, 8 = 0.098, 
E52 = 50.5, Ay = 0.064 mag, and reduced chi-squared value 42.1/26 = 1.62. All light curves show 
the total of synchrotron and inverse Compton emission; optical and near-infrared data have been 
converted to J-band flux densities for clarity in plotting. For the best-fit model only, we plot the 
contribution to the X-ray flux from inverse Compton emission separately, as the red dash-dotted 
line. Bottom: The chi-squared contribution from each data point. Positive values indicate that 
the best-fit model underestimates the flux, and negative values indicate that it overestimates the 
flux. 
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Fig. 5. — Posterior distribution of all parameters in model (B), in which the flares at t m 470 s are 
considered to be emission from the reverse shock regions: p, the power-law index of the electron 
energy distribution; e e , the electron equipartition parameter; €bj, the magnetic field equipartition 
parameter in the forward shock; n, the circumstellar density in units of cm -3 ; 9, the opening half- 
angle before jet break; E52, the isotropic-equivalent kinetic energy in units of 10 52 ergs; Ay, the dust 
extinction of the host galaxy; Tq, the initial Lorentz factor of the outflow; and €B, r , the magnetic 
field equipartition parameter in the reverse shock. The shaded blue region delimits the 1-cr range 
(68.2%), and the region included within the vertical blue lines corresponds to 90%-confidence level. 
The green color indicates the region for j m < 2.1 (the red color corresponds to 1-cr confidence 
range), which is around 4.7% of the total trials; ; its height has been magnified by a factor of five 
for clarity. 
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Fig. 6. — Top: Theoretical light curves (solid lines) corresponding to the best-fit parameters in 
model (B), in which the flares at t ~ 470 s are considered to be emission from the reverse shock 
regions. Best-fit model parameters are: p = 2.243, e e = 0.0084, ebj = 5.7 x 10 -3 , n = 212. 4cm -3 , 
6 = 0.126, E 52 = 147, Ay = 3.18 x 10" 2 mag, T = 183.6, and e B , r = 0.50. The dotted lines 
indicate the separate flux contributions from the reverse and forward shocks (reverse shock emission 
is distinguished by its fast decay at early times). The solid lines indicate the total model flux, 
with the red dash-dotted line showing the contribution to the X-ray flux from inverse Compton 
emission, which is relatively unimportant compared to the synchrotron component. All optical and 
near-infrared data have been converted to J-band flux densities for clarity in plotting; underlying 
model calculations use the various observed frequencies directly. Radio, /-band, J-band, H-b&nd, 
X-ray, and BAT data are plotted in blue, orange, green, yellow, red, and magenta, respectively, as 
indicated in the plot legend. Bottom: The contribution to the total (best-fit) model chi-squared 
for each data point. Positive values indicate that the model underestimates the flux, and negative 
values indicate that it overestimates the flux. 
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Fig. 7. — Joint confidence regions (68%, 90%, and 99%, respectively) for three chosen parameters 
from our model fits: the jet opening half-angle 9, the isotropic-equivalent kinetic energy E52, and 
the circumburst density n. Top: Joint confidence regions for model (A), which excludes data 
related to the flaring activity at t ~ 470 s. Bottom: Joint confidence regions for model (B), in 
which the flares at t ~ 470 s are considered to be emission from the reverse shock regions. With the 
additional constraints available in model (B), the separate effects of blastwave kinetic energy and 
circumburst density can be distinguished, and the final constraints on the beaming of the burst 
and its total kinetic energy are significantly better-defined. 



